*
/This is SAS syntax for nonlinear Model 2: bifactor with 1 group factor             **
/The derivation of the nonlinear SEM reliability and the original SAS codes         **
/are presented in Green & Yang (2011b). Polychoric correlation, loadings,           **
/and thresholds are obtained by fitting the model in Mplus with WLSMV method.       **
*************************************************************************************;


proc iml;
RESET fuzz;
THRESH={
 -1.307  -0.677   0.090 ,
 -1.144  -0.604  -0.073 ,
 -1.112  -0.688   0.039 ,
 -1.823  -1.307   0.021 ,
 -1.050  -0.327   0.176 ,
 -1.572  -0.861  -0.309 ,
 -1.439  -0.656   0.345 
};

LOAD={
0.532 .793,
0.460 .547,
0.557 .710,
0.800 0,
0.686 0,
0.874 0,
0.636 0
};

FACCOR={1 0,0 1};

POLY={
1       0.679    0.860    0.421    0.297   0.471   0.410,
0.679  1         0.645    0.317    0.323   0.485   0.212,
0.860  0.645    1         0.482    0.330   0.422   0.426,
0.421  0.317    0.482    1 	   0.516   0.684   0.565,
0.297  0.323    0.330    0.516    1        0.662   0.403,
0.471  0.485    0.422    0.684    0.662   1        0.455,
0.410  0.212    0.426    0.565    0.403   0.455   1
};

NTHRESH=Ncol(thresh);
NCAT=NTHRESH+1;
NITEM=Nrow(LOAD);
NFACT=Ncol(LOAD);

POLYR=LOAD*FACCOR*T(LOAD);print polyr;run;
do j=1 to NITEM;
POLY[j,j]=1;
end;
DIFFPOLY=POLY-POLYR;

Print NTHRESH[label="Number of Thresholds"], NITEM[label="Number of items"], 
NCAT[label="Number of response categories"], NFACT[label="Number of factors"], 
THRESH[label="Response Thresholds"],LOAD[label="Factor Loadings"],
FACCOR[label="Factor Correlation Matrix"],
POLY[label="Polychoric Correlation Matrix among Continuous Items"] ; 
print "The matrix below is the difference between polychoric correlation matrix generated
 by factors and inputted polychoric correlation matrix. Nonzero values should 
 represent the estimated correlated errors, as specified by the user, or an error in inputted data.";
print DIFFPOLY[label=" "];


/* Computing numerator and denominaot of Equation 21 */ 
sumnum=0;
addden=0;
do j=1 to NITEM;
do jp=1 to NITEM;
sumprobn2=0;
addprobn2=0;
do c=1 to NTHRESH;
do cp=1 to NTHRESH;
sumrvstar=0;
do k=1 to NFACT;
do kp=1 to NFACT;
sumrvstar=sumrvstar+LOAD[j,k]*LOAD[jp,kp]*FACCOR[k,kp];
end;
end;
sumprobn2=sumprobn2+probbnrm(THRESH[j,c],THRESH[jp,cp],sumrvstar);
addprobn2=addprobn2+probbnrm(THRESH[j,c],THRESH[jp,cp],POLY[j,jp]);
end;
end;
sumprobn1=0;
sumprobn1p=0;
do cc=1 to NTHRESH;
sumprobn1=sumprobn1+CDF('NORMAL',THRESH[j,cc]);
sumprobn1p=sumprobn1p+CDF('NORMAL',THRESH[jp,cc]);
end;
sumnum=sumnum+(sumprobn2-sumprobn1*sumprobn1p);
addden=addden+(addprobn2-sumprobn1*sumprobn1p);
end;
end;
reliab=sumnum/addden;
print sumnum[label="Numerator of Equation 21"], 
      addden[label="Denominator of Equation 21"],
      reliab[label="Nonlinear SEM Reliability Coefficient"];



*
/The following SAS codes are written based on the general SEM method for computing reliability **
/in Raykov and Shrout (2002). Ftrue1 to Ftrue3 present true score variance due to F1 to F3.    **
/Ftrue and Ftotal is the scale true score variance and scale score variance, respectively.     **                                             
 

/This is SAS code for linear Model 2: bifactor with 1 group factor                  **
/This model yielded variances of true and total of 19.37 and 22.53, respectively.   **	
/Reliability is computed by 19.37/22.53=.860                                        **
**************************************************************************************;

data a; set demodata;run;
proc calis DATA=A method=max cov maxiter=2000; 
lineqs 
x1=a1 F1 +b1 f2  +e1,
x2=a2 F1 +b2 f2  +e2,
x3=a3 F1 +b3 f2  +e3,
x4=a4 F1 +e4,
x5=a5 F1 +e5,
x6=a6 F1+e6,
x7=a7 F1 +e7,

Ftrue1=g1 F1 +d1,  
Ftrue2=g2 F2 +d2,
Ftrue=Ftrue1+Ftrue2 +d3, 
Ftotal=x1+x2+x3+x4+x5+x6+x7 +d4; 

std 
e1-e7 =ve1-ve7,
f1=1,
f2=1,
d1=0,
d2=0,
d3=0,
d4=0;

cov
f1 f2=0;

lincon
g1-a1-a2-a3-a4-a5-a6-a7=0,
g2-b1-b2-b3=0;
run;


*
/This is SAS syntax for linear Model 3: bifactor with 2 group factors               **
/This model yielded variances of true and total of 19.72 and 22.51, respectively.   **	
/Reliability is computed by 19.72/22.51=.876                                        **
**************************************************************************************;
data a; set demodata;run;
proc calis DATA=A method=max cov maxiter=2000; 
lineqs 
x1=a1 F1 +b1 f2        +e1,
x2=a2 F1 +b2 f2 +c1 f3 +e2,
x3=a3 F1 +b3 f2        +e3,
x4=a4 F1               +e4,
x5=a5 F1        +c2 f3 +e5,
x6=a6 F1        +c3 f3 +e6,
x7=a7 F1               +e7,

Ftrue1=g1 F1 +d1,  
Ftrue2=g2 F2 +d2,
Ftrue3=g3 F3 +d3,
Ftrue=Ftrue1+Ftrue2+Ftrue3 +d4, 
Ftotal=x1+x2+x3+x4+x5+x6+x7 +d5; 

std 
e1-e7 =ve1-ve7,
f1=1,
f2=1,
f3=1,
d1=0,
d2=0,
d3=0,
d4=0,
d5=0;

cov
f1 f2=0,
f1 f3=0,
f2 f3=0;

lincon
g1-a1-a2-a3-a4-a5-a6-a7=0,
g2-b1-b2-b3=0,
g3-c1-c2-c3=0;
run;


*
/This is SAS syntax for linear Model 4: bifactor with 1 group factor               **
/and correlated errors among item 2,5,and 6                                        **
/This model yielded variances of true and total of 18.12 and 22.51, respectively.  **	
/Reliability is computed by 18.12/22.51=.805                                       **
*************************************************************************************;
data a; set demodata;run;
proc calis DATA=A method=max cov maxiter=2000; 
lineqs 
x1=a1 F1 +b1 f2  +e1,
x2=a2 F1 +b2 f2  +e2,
x3=a3 F1 +b3 f2  +e3,
x4=a4 F1         +e4,
x5=a5 F1         +e5,
x6=a6 F1         +e6,
x7=a7 F1         +e7,

Ftrue1=g1 F1 +d1,  
Ftrue2=g2 F2 +d2,
Ftrue=Ftrue1+Ftrue2 +d3, 
Ftotal=x1+x2+x3+x4+x5+x6+x7 +d4; 

std 
e1-e7 =ve1-ve7,
f1=1,
f2=1,
d1=0,
d2=0,
d3=0,
d4=0;

cov
f1 f2=0,
e2 e5=cov1,
e2 e6=cov2,
e5 e6=cov3;

lincon
g1-a1-a2-a3-a4-a5-a6-a7=0,
g2-b1-b2-b3=0;
run;